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ABSTRACT 


The "Fast Fourier Transform" (suggested by Drs. 
J. W. Cooley and J. W. Tukey) is presented with special 
application to solving the interferogram integral obtained 
in Fourier Spectroscopy. Computational timing is tabu- 
lated and an explanation of a computer routine using this 
method to a binary base is presented. 
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THE 

FAST FOURIER TRANSFORM 
IN 

FOURIER SPECTROSCOPY 


INTRODUCTION 

The problem of digitally reducing interferogram data obtained from an in- 
terferometer spectrometer is two fold; (1) numerical quadrature of the inter- 
ferogram integral to obtain the spectral magnitudes, and (2) computation of the 
phase of the spectral magnitudes to determine the direction of radiation flow. 

A third requirement can be placed on the computation if a great number of 
spectrums are to be analysed; that is fast calculation on a high speed digital 
computer. This is the problem with which this paper is primarily concerned. 

Many methods of fast calculations of the interferogram integral have been 
suggested; however, the "Fast Fourier Transform" (Cooley, J. W. and Tukey, 

J. W., 1965) method offers many advantages in speed and accuracy which others 
do not have. 

The application of this method to discrete interferogram data is discussed 
and the use of an existing computer program written by Dr. J. W. Cooley in 
FORTRAN IV using the "Fast Fourier Transform" to compute a Fourier trans- 
form or series is explained in Appendicies I and II. Only a small amount of in- 
terferometer theory is discussed so as to keep the general theme of the report. 

It should be mentioned that the Fast Fourier Transform has been imple- 
mented with much success into the Infrared Interferometer Spectrometer (IRIS) 
Experiment Data Reduction Program as well as the theoretical simulating the 
IRIS experiment. Computational timings have been tabulated using these pro- 
grams and are given in a following section. The reader should also be aware 
that now with the computational speeds obtained using the "Fast Fourier Trans- 
form," problems in Fourier Spectroscopy which previously were overly time 
o nisuming on the computer are now realistic for computer solution. These in- 
clude taking convolutions for truncated Fourier Integrals, shown by Dr. J. W. 
Cooley, co-author of the Fast Fourier Transform (in an unpublished report) and 
correction of asymmetric interferograms (Forman, M. L., 1965). 


1 


The author is greatly indebted to Drs. B. Conrath and R. R. Hanel for many 
hours of helpful discussions on interferometer theory and Fourier analysis, and 
to Dr. J. \V. Cooley for the use of his computer program. 


DESCRIPTION OF THE PROBLEM 

The instrument used to obtain interferogram data in the IRIS experiment is 
essentially a two beam Michelson interferometer, sketched in Figure 1. 



DETECTOR 


Figure 1— Sketch of Michelson Interferometer 


The incoming heterochromatic radiation is split at the beam splitter, A, 
into two waves of equal amplitude; one directed toward a moveable mirror at C, 
and the other toward a fixed mirror at B. They are reflected back to the beam 
splitter, recombined, and directed to the detector at D. The recombined signal 
received at the detector is the interference pattern of the two beams called the 
interferogram and, ideally, can be defined as a function of path difference traveled 
by the two beams by 


1 ( 8 ) = 



The path difference, S, is defined with the aid of Figure 1 by 


( 1 ) 


S = 2 (CA - BA) 


2 



The multiplying factor of two times the distance (CA-BA) arises from the fact 
that the radiation travels to the mirror and back again to transverse the added 
distance (8 can be given as a function of time and mirror velocity; however, we 
will use the path difference relationship here). 

is the spectral distribution of the incoming radiation and is written here 
in terms of wave number, v . It has been extended to include negative wave 
numbers by the definition 


B 


B. 


Clearly, B y is real, and B v and 1(8) are transform pairs; thus one is able 
to obtain the spectrum B„ by taking the Fourier transform of equation (1) which 
is 

/• CD 

B v = 1(8) e~ i2 ™ s d8. (2) 

•'-00 

If 1(8) were symmetric, which in the theoretical case is true, equation (2) 
reduces to 

B v = 2 f 1(8) cos (277 v 8) dS 

vo 

However, in actual interferometer use a phase shift, <p , is introduced by 
the instrument components to the incoming radiation and equation (1) should be 
re-written as 


I(S) = 



B v e ltPv e i27rvS dv 


(3) 


1(8) in this case will be an asymmetric function about zero path difference, 
8 . Upon taking the Fourier transform of equation (3), one obtains 


i® f ” 

B v e v =\ l( 8 

» -m 


1(8) e ~ i27TvS dS. 


(4) 
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Equation (4) is of the form 


B v e icpv = Cv + i Sv 


and one merely takes the absolute value of the right hand side to obtain the 
spectral magnitudes. 

The phase shift cp v can, of course, be computed from 



The phase spectrum can be a useful tool in interpreting the amplitude 
spectrum. The direction of the net energy transfer between the detector and 
the target (determining whether the detector is warmer or colder than the 
target) changes direction from one spectral region to another and is indicated 
by an abrupt phase shift of 180°. 

Therefore, the problem of digitally reducing interferogram data requires 
solution of the interferogram integral in equation (4) and computation of the 
spectral amplitude phases by equation (5). 


NUMERICAL QUADRATURE OF THE INTERFEROGRAM INTEGRAL 


^ The numerical calculation of the integral in equation (4) can be done by 

various quadrature methods. The Gaussian quadrature is perhaps the most 
accurate, but this method requires unequal spacing in the sampled data. 

It is suggested that in practical use, when the truncated range (- S lt S 2 ) of the 
integral (see equation (6)) is large enough, many advantages are offered by use 
of the trapezoidal rule. 


One advantage in accuracy is easily seen by looking at Eulers summation 
formula: (Scarborough, 1955) 


r 


f(x) dx = ‘h 


f ( x o> 


f(*n) 


7 20 


2 + f(X l ) + • • • + f ( X n— 1 ) + 2 

[f v (b) - f v (a)] 


-f 2 


+ _ — b 5 


30240 
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This states that if the odd derivatives at the endpoints of an integration are 
small or equal, the trapezoidal rule is an excellent approximation to the integral. 

In many cases, the derivatives at the endpoints of the interferogram are 
equal, and in the worst cases observed the derivatives were always of the 
same order. In these cases, the order of the error would be h , the separation 
distance (AS in the interferogram integral), and for the interferogram integral, 
this is small. 

Therefore, the trapezoidal rule appears to be a satisfactory quadrature 
method and, as one will see, is a suitable form for the "Fast Fourier Trans- 
form" to obtain speed of calculation. 

Rewriting equation (4) as a truncated inverse over the interval (-5 , 8 2 ), 
and, for convenience, setting =B v e iCf> v , we have 


£ 


K = I We 
. 8 , 


- i 27rvS 


( 6 ) 


Since the interferogram values sampled by the instrument are an average 
over an interval, the trapezoidal rule takes the form 


A v = B y e = AS 


Z] i<s i )e “ 


- i27TvS j 


(7) 


J --m 


where the subscript, j , ranges over N sampled intervals with N defined as 


N = m + n + 1 


The point 8 is of course the sampled interval which normally will not 
correspond exactly to the point of zero path difference. Correction methods 
have been studied for correction of asymmetric interferograms of this type 
(Forman, M. L.); however, this paper will not concern itself with these methods. 

Solution of equation (7) to obtain the spectral magnitudes and phase angles 
is easily accomplished using the "Fast Fourier Transform," and, as one will 
see, yields a very high degree of accuracy for computer solution. 
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THE "FAST FOURIER TRANSFORM" 


1. Description of the Method 

The form of the Fourier series equation required for the "Fast Fourier 
Transform" is 


I(j) 



k = 0 


• 277 i k 
B(k)e 1 - N 


( 8 ) 


j = 0, 1, 2. . . N-l 


and its transform is 


N - 1 . 2-rrj k 

B(k) = 1. ^ I ( 8 j) (9) 

j =0 

Normal solution of equation (8) using a decimal based summation requires 
an order of N 2 operations where an operation is defined as one multiplication 
and addition. However, the number of operations using a different base, say r , 
is of the order 8N log r N . A simple example will illustrate the saving. 

If N is composite such that N =r x • r 2 , and j and k are written 

j = Tj + j q , }q — 0, 1, 2, • • • Tj — 1; j j -0,1, • • • r 2 1 


and 


k = r 2 kj + k Q> k Q = 0, 1, . . . r 2 - 1; k x = 0, 1, ... r j - 1 
then equation (8) can be written 


6 


i(j) = 7 -1 


' *Tr j (f 2 k l + V 


L \ L i 

k o k l 


B(k) e 




B(k) e 


t ^ jr k iil jkn 

N J 2 1 N 0 


k 0 k l 


But 


2 7T . . 2 77 2 7 7 

1 N ' '2 k l = e A — k l r 2 < r l j l + j 0) k l i 0 


Therefore, the inner sum over k x , depends only on j Q and k Q , which can be 
written as 


B i (io' k o) = B ( r 2 k i + k 0 ) e 


. 2 77 

7T k i J o r 2 


and equation (8) now becomes 

r ( r X i i + j 0 ) = ^ B ! (j 0 - k o) e 


277 

1 — j k o 

N 0 


The total number of operations for this case is now N(r x + r 2 ) as opposed to N 2 before. 

It has been shown (Colley, J. W., Tukey, J. VV., 1965) that the most saving is 
obtained when N is written as some number raised to an integer power, or 


N = r m 
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and, further, if r = e, one uses the least operations in solving equation (8). How- 
ever, if r = 2 , the saving is approximately the same, and from a digital com- 
puter standpoint, certain advantages in programming are obtained. 


With this in mind, if one writes 


N = 2 m 


and 


i = 2 m ~ 1 j , + 2 m " 2 j „ + 

J J m - 1 J m-2 


• + i 0 


k 


k m- 1 + 2 m- i 


+ 


k 


0 


where each and k^ , -i = o, • • • m-1 take on the values 0, 1, then equation (8) 
takes the form 


1(2 


m-l: + 2 m “ 2 i + 

Jm.l + 2 J m-2 T 


+ i 0 ) = B .a„ + 2 j, + ...+ 2-> j..,) 


where 


B i(j 0 , 2 m - 2 k m _ 2 + . . . +k 0 ) = B ( 2 ' m ' lk -i 

k m - 1 


+k 0> e 


i^j 0 k 1 2“- 1 

jsj * 0 m — 1 


8 


and 


4.,, k 0 + ... 

]2 H.i (io + ••• + 2 4 - 2 U-r 

k m - / 6 


277 - 

1 IT * A- i 
e 


2 


-t- l 


+ 


+ in) 


,a 


A 


/ C = 1, 2, . . . , m 


The routine listed in Appendix I uses the base two and is written such that 
either the Fourier series coefficients or their transform can be computed. That 
is, one is able to compute either equation (8) or equation (9). 


2. Application to the Interferogram Integral. 

We saw that the spectral magnitudes and phase angles can be obtained from 
equation (7), 


k v - By e iCpv = AS 1(8.) e -i2lTvS i. 


j = -m 


If one sets, 


v k = kAv, k = 0, 1, . . . N - 1 


and 


8 = j AS, j = 0, 1, . . . N - 1 
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where v Q = S 0 = 0 , we have 


Ay, - AS 

k 


m 

£ ^ 

j =- n 


g -i 2 rr j k A vAS. 


Further, if- one makes the substitutions 


( 10 ) 


j = i - n 


ASAy = N" 1 = (Number of sampled intervals)" l . 


one can write the summation in the form suitable for the "Fast Fourier Trans- 
form." That is, 


Ay = AS 

k 


2rrkn jsf 1 



3 = 0 


2 TT 

- i j k 

I(S.)e N 


( 11 ) 


where 


N = m + n + 1. 


If the number of sampled intervals, N , does not equal a power of two and 
the interferogram has zero mean (no DC component), one is able to add values 
equal to zero to either or both ends of the interferogram without loss of ac- 
curacy and compute spectral magnitudes corresponding to v =kAy. Here, 


Ay 


1 

AS • N' 


where 


N' = N + l 


such that N' is equal to two raised to some integral power. 
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Figure 2-Interferometer extended with zeros on each end with 
S Q centered in area of N* intervals. 


Further, if the interferogram is positioned in the area of N' intervals as 
shown in Figure 2, such that S Q is a distance N'/2 from the beginning, a par- 
ticularly simple case arises. Keeping in mind that A v = > equation (11) 

takes the following form for this case, 

N ” 1 2 7T 

- i — j k 

Av k = AS e i7Tk ) 1 ( 8 .) e N (1 

i = o 

and the spectral magnitudes are computed as before from 


Bv k = |A^ k 


(13) 


Re-writing equation (12) as 


Bz^, e 


i cpy 


k _ 


k 7 T 


( Cu k ~ i S vk> 


or 


Bv k g 1 (^k " k7r > = C u k 


i S 


vk 


the phase angle can then be computed from 


qpv. = tan 


-‘l- 

\Cv, , 


+ k 77 


(14) 
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or 


<pv k 


tan " 1 


(-l/Si-i, 

(-l) k Cv k 


For the case when the interval § 0 is not positioned in the center of the array, 
one still can use equation (12) to compute the spectral magnitudes, however, one 
must use an equation similar to equation (11) to compute the phase angles. That 
is, 





27/kn 

N 


(15) 


where, n , is the number of intervals from the beginning of the array to S 0 . 

In using the routine in Appendix I, one must be careful applying it to com- 
pute the Interferogram integral. From experience it proved easier and faster 
to compute the Fourier series coefficients rather than the Fourier transform. 
This means that for the real input data, one obtains the complex conjugate of the 
vector B?/ k e lCpl/k as output from the "Fast Fourier Transform." 

Clearly, this makes no difference in computing the spectral magnitudes 
other than a factor of two, however, the phase angles must now be computed 
from 





(16) 


or 


q>v k = tan 


‘(-l) k + 1 Si 


(- 1 ) 


k + 1 


Cv, 


k J 


(17) 


or, in the case where 8 is not positioned in the center of the interferogram, 


cp^ = 


27ikn 

N 



(18) 
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3. Summary of the Application. 

The method of computing the spectral magnitudes and phase angles can be 
summarized as follows: 

1. Read the interferogram data into the "Fast Fourier Transform" (FFT) 
and compute a Fourier series yielding Cv k and Su k as output. 


FFT (series) = Cv k + i S^ k 


2. Compute the spectral magnitudes from 


= 2 • AS • |Ci/ fc + i S^ k 


where 


and 


v k = k Ay, k = 0, 1, . . . N - 1 


Av 


1 

AS -N' 


3. Compute the phase angles from either equation (17) or (18) depending 
upon the position of the interferogram data in the array of sizeN' = 2 m . 

In using the routine in Appendix I, one automatically obtains the same 
number of points of output as he has read in as input. 

4. Computational Timing 

Figure 3 shows the computational times obtained using the FORTRAN IV 
version of the "Fast Fourier Transform" routine (HARM) listed in Appendix I. 
The line indicating the times obtained on the IBM 360/65 system were obtained 
with the HARM subroutine, compiled with the optimization option. These times 
were obtained by Mr. Guy Marcot, Laboratory for Space Sciences, GSFC. A 
comparison with conventional times obtained with the algorithm by Goertzel 
indicates the saving, especially with a large number of data points. 
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Figure 3— Computation Times using the FORTRAN IV version of HARM on 
the IBM 7094 and 360/65 Computers vs. the Goertzel algorithm. 

Another version of the HARM subroutine, written in FORTRAN II and FAP, 
which is used in the theoretical simulation of the IRIS experiment was timed for 
2 1 A and 2 12 points on the IBM 7094 computer. These results are tabulated 
below. 


No. of 

Time for Fourier 

Points 

series calculation 

2 1 1 

1.3 sec 

2 12 

2.6 sec 


The times tabulated above include the time obtained to compute the sine and 
cosine transform and do not include taking the absolute value or obtain the phase 
angle. This has been timed using the theoretical IRIS simulation. To compute 
722 spectral magnitudes and phase angles out of the 2 11 or 2 12 available took 
approximately 3/10 seconds. 
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APPENDIX I 


Sample Program Using the "Fast Fourier Transform" Routine "HARM" 

The following program has been written exemplifying the use of the "HARM" 
subroutine to compute the transform of a set of real data and take the inverse 
of the transform to compare with the original data. 

"HARM" is written in FORTRAN IV and is set up to accept complex input as 
well as to perform a three-dimensional sum in directions M(l), M(2), M(3). 

The return from "HARM” has the complex vector computed and stored in the 
array A, where 


A(l) + i A(2) = Cv Q + i Sv Q 
A(3) + i A(4) = + i Sz>j 


and so on. The array M and option IFS must be set prior to entry, and if one 
wishes to use the routine in one -dimens ion, as in the sample program, one sets 


M( 1) = log 2 N 
M(2) = 0 
M(3) = 0 


The call to the "HARM" subroutine is 

CALL HARM (A, M, INV, S, IFS, IFERR) 

A = Array of complex input, where the real and imaginary parts are 
in consecutive storage locations. A must be dimensioned 2N. 

M = Array containing the logrithm to the base 2 of the length of the 
summations in directions M(l), M(2), M(3). M is dimensioned 3. 
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INV = Array computed in "HARM" for bit inverting. INV is dimensioned 
N/8. 

S = Array computed in "HARM" containing the array of sines. S is 
dimensioned N/8. 

IFS = 0 compute INV and S tables 

+1 compute INV and S tables and do Fourier series 
-1 compute INV and S tables and do Fourier transform 
+2 do Fourier series only 
-2 do Fourier transform only 

IFERR = Error option. (See Appendix II) 

Further amplification of the use of "HARM" can be obtained from the 
program write-up on file in the Laboratory for Atmospheric and Biological 
Sciences. 


18 



c. TEST OF FAST FOURIER TRANSFORM ROUT I NE ( HARM ) • 

C THIS PROGRAM PERFORMS A TRANSFORM ON A SET OF DISCRETE REAL DATA » 

C THEN TAKES THE INVERSE Of THE TRANSFORM TO RECOVER THE ORIGINAL 

' C DATA. THL INPUT DATA. REPRESENTS A BLACK BODY TEMPERATURE SPECTRU! 

C COMPUTED IN THE 'IRIS' EXPERIMENT. 

DIMENSION A( 16384) »M(3 ) > I NV ( 2048 ) »S( 2048) 

C SET UP DIMENSIONS 

MCI) =6 
M ( 2 ) =0 
M ( 3 ) =0 
N 1 = 2**M(1) 

N 2-2 - r ' M ( 2 ) 

N3 = 2**M( 3 ) 

NT0T=N1*N2*N3 
DO 10 I=1,NT0T 
10 AID =0. 

READ (2,15) ( A ( 2 * I - 1 ) » I = 1,NT0T) 

15 FORMAT (6 £12. 5) 

16 WRITE (3,17) ( A ( 2 * J - 1 ) , J = 1 , N T OT ) 

17 FORMAT (26H10RIGINAL DATA- REAL PART //(5F13.5)) 

WRITE (3,18) (A(2*J), J=1,NTGT) 

18 FORMAT (26HOORIG INAL DATA- I MAG PART //(5F13.5)) 

C EVALUATE FOURIER TRANSFORM 

C IF IFS = -1, SETUP INV AND S TABLES AND AND COMPUTE FOURIER SERIES 

IFS = -1 

CALL HARM ( A , M , I N V , S , I FS , I FERR ) 

25 WRITE ( 3 , A 2 ) ( A ( 2*J- 1 ) » J= 1 »NTOT ) 

42 FORMAT ( 23H1REAL PART OF TRANSFORM/ / ( 5F13.5 ) ) 

WRITE (3,44) ( A ( 2*J ) ,J = 1 ,NTOT ) 

44 FORMAT ( 15HG IMAGINARY PART / / ( 5 F 1 3 . 5 ) ) 

C TAKE INVFRSE 

C IF I FS=+2 , COMPUTE FOURIER SERIES. 

IFS = +2 

CALL H A R M ( A , M , INV,S, IFS, IFFRR ) 

WRITE (3,70) ( A ( 2*J-1 ) , J=1 ,NTOT ) 

70 FORMAT ( 60H1REAL PART OF ORIGINAL INPUT OBTAINED BY INVERTING TRANS 
1F0RM// ( 5F13 .5 ) ) 

WRITE (3,71) ( A ( 2*J ) ,J = 1 ,NT0T ) 

71 FORMAT ( 15H0IMAGINARY PART / / ( 5 F 1 3 . 5 ) ) 

PAUSE 77777 

END 
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ORI GTNAL DATA- RFAt. PART 


23C.0<;00C 229* 99000 229. 99C00 230*01000 

ZZC+JIBSAL 23C *JL8£.Q3l..- — .230*. 0 200 C — . 

23C.030CO 230. 10 000 229.96000 229.95GCC 

?;MV100C<> 229* 94000 - 229* 98C00 23C. 03CO0 

230*01000 229. R4000 229. 99C90 230*07000 

.?£9*«?0£C - -2-29*-85£n0 230* 00000 229. 90000 

23C.C50PC 229. 88000 229.84000 229. 4700 C- 

p.g.Q^-tWuTCr; 22 9* 5 800.0 2 29 * 8.90010 23Q*Q.4QQC. 

229. 8 5000 229. 74000 229.61000 229* 21300 

229« 52000- 2-29*-2 7000 — 229* 630,0 0 - 229*69000 

229* 10000 229, 71000 229.63000 229*24000 

- --- ? 29* 42000- — 22 9* 5-7000 - - 229* 5900 0 229* 950 QC 

230. 00000 22 9.87000 229.67000 229. 64CC0 


ORIGINAL DATA- IMAG PART 


0. 0* Oc 9. 

— c* — o* o* - — o. 

0. 0. 0. 0. 

a. a. o* a* 

C. 0. 0, 0. 

- — 0 # 0* 0* — - . o 

0. O. -0.00000 -0*00000 

-C. ronro -0*00000 =0. 00000 tO.GOOQC .. 

-c.ooooo -o.ccooo -c.ooooo -o.oocoo 

= 0*00000 t Q. o op.oc -c.ooono^- =0 *.aiooji_ 

-Cm 00000 cooor -c.ooooo -c.cocoo 

-0* 000 CO -XU 00000 -0*00003 -0*00 a 03 

-r.roorn -o.cgooo - 0.00000 -o.ooocc 


230.09003 

22 . 9 .- 903-03 

229.99003 
230.C8000 
229. 98003 
229*74003 
229. 73003 

..-.23-3*4 60-u 0 

229.80003 
229.09000 
229.05300 
229*88003 


0 » 

.... 3, 

3. 

~ 3. 

Oft 

-0.00000 
-3.C0003 ... 
-O.COOOO 

--o...oaoOL'3 

-0.00003 
- 0*00000 
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R F At PART CF TRANSFORM 


?29„ 81062 Or PQ 5 26 0. 04921 0.03292 -0.01639 

— C* 0C.3CJ -£~£33m riWQL32,6. -C. 00515 -Q. 01372 

C.02211 0.02160 -0.00379 -0.00196 -0.00677 

— 0. 0135 9 0*0 2375 0.01 425 -0. 02524 0 .0131 1 

0.02518 -Or. 00713 0.00676 0.01149 0.00546 

- C Q 00.945 .0.-015-3.9 -0.0185 4 0.00375 0. CO 94 5 

C.C09C7 0.00188 C. 01375 0.00188 0.00907 

Ca.C.2945 .. -Jlt.0JL8.5_4_-. 0.015 39 ... -0, PC 9 45 

Op 005 46 0.01149 0.00676 -0.00713 0.02518 

— .-..r*Q13. 1J -O r 075 24 0 101425 0. 02375 0.01359 

-0*00677 -0.00196 -0.00379 0.02160 0.02211 

-0.0 1372 -o,. 00515 -0.01 326 -O.C 3 302 0.00301 

-0.01639 0.03292 0.04921 0.00526 


IMAGINARY PART 


0. 

-0. 1 3424 

-0» 03590 

0.02090 

0.00464 

Op 00 215 

Oft 014 5.9 

-Or. 02 89 4 

0.00230 

-0.03026 

-0*01 345 

0. 0 0272 

0« 0044 9 

-0.00504 

-0.01928 

-n^03Aaa 

... Oa. c 3 250 

. 0. 0128 5 

Oo 00 1 26 

-0.04073 

0.00637 

0. 0 05 57 

-0. 00 82 7 

0 . 0046C 

0.00698 

Co 0.0 565 

-0«i 00.509 . ...... 

0,00575 

-0,01536 

0.00598 

-0.00021 

-0. 00435 

0. 

0« 00435 

O.C 00 2 1 

-0.00598 

,0.01536 

-0. CO 57 5 

0.00509 

-0.00565 

-0.00698 

-Or 00460 

0.00827 

-0.00557 

-0.00637 

CVOA ail 

-.<3a.Q11.26. 

-Do 012.8 5 

-0,03250 

O.C 3480 

0. 01 928 

0.00504 

-0. 00449 

-0,00272 

0.01345 

0.03026 

— o. 0 0230 

0 , 0 2 8 9 4 

-Oft 01459 

-0.00215 

— 0* 03464 

-0. 0 20 90 

0. 03590 

0. 13424 
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REAL PART CF GRIGTNAL INPUT 

OBTAINED BY 

INVERTING TRANSFORM 

230.0 8966 

2 2 9, 9 8999 

229. 98999 

230,00999 

230,08999 



J|Bgp?f2PYiWZ*T*T*HH 

2^. Til 999 

229. 89999 

23o„ r?666 

230,09999 

229, 95999 

229,94999 • 

229,58999 

- 9-9^9 

2 29, 93999 

22 9* 97-9-9.9. 

-230,02999 

.230,079.99 

?3<% C0666 

229, 83999 

229. 98999 

230,06999 

2 29,97999 

- - 2 2 6. -86 9 59 

229, 8499.9 

.229. 9999 9 

229^.89.9.9.9 

P79* 72999 

230,04669 

229, 88000 

229. 83999 

229, 46999 

229.72999 

-^r 91999 

220^5-7^99. 

22 9* Rft.99 9 

230. 03969 

? 2H- 0 5999 

229, 84999 

229, 73999 

22 9, 6099 9 

229, 21000 

229. 79999 

2 29^ 519 -9 9 

-_J?_29^26 999 - 

22S*~62999 

.... 229 * 68999 

? ? 9 * 0 899 9 

229. 09699 

229, 70999 

229. 62999 

229,23999 

229.04999 

- -- .. 2?9« 44 999 

- 2-? Q<f-56-999 

. 22 9. 5399.9 

229.94999 ._ 

2P9 » 8 7999 

??9„ 99999 

229, 86999 

229. 66999 

229,64000 


IMAGINARY PART 





o. 

-Oo ^OOOC 

-0,00003 

-0,00000 

0. 

— . c*.nnn.oo . . . . 

._o...Goam 

.. -0.0303 3 

o. 

3.00000 

c* nnoco 

-o, on ooo 

0. 

-0,00000 

-3,00000 

. 

. r>. 



0.C3000 

0 . 

o,ooono 

0« 

0.00000 

0. 

- e. ooaea— .. 

— -jo, ono on. 

f% 

.. ....... Gu... 

0.03000 

r, on 9 ro 

o, ccoon 

o * 

0,03000 

O.OCOOO 

- ocnncn ... 

r lf . 

-0.09000 

-0.00000 

3.00003 

r. 

G, OOGfiO 

-0. 00000 

O.GOOCC 

D* 

_ r. onnnn 

o, oo - or. 

BSPCnSTilrolK 

0. 

-3,00000 

-C. 00000 

-C. 00030 

0. 

-3, GOOOO 

0.C00C3 

-g. no nor 

0* 

-n, nnono 

0.03003 

0. 

r. 

-r> oaooo 

-C. 00000 

-0,00000 
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APPENDIX H 

"FAST FOURIER TRANSFORM" - "HARM" 


C HARM DISCRETE FOURIER TRANSFORM. BASIC FORTRAN IV, IBSYS. 

simeruTi ne harmia.m, inv , s, ifs,. IE err ) _ 

DIMENSION A C 1 > , INV(l) ,S(1 ) , N( 3) ,M( 3 »,NP( 3) ,W(2 ) ,W2< 2 ) , W3(2 ) 
EOm-VAi F N CE . IN i4t I U , ( N2. N(21 1 , IN3 , N (3 ) I _ 

— C, -INPUT, ■e-AHAMETERS.-ia.-aF - SEI-BX- USES B££3££- ENTERING HARM- 

£ A IS A 3—0 1 MF NS IONA I ARRAY OF COMPLEX COEFE ICIENTS, 

0 OF OIMFNSION(Nm,N<2),N(3n* 

- C THE A*S ARF STORED WITH REAL PART OF A(I1. 12. 13) IN THE LOCA TION 

0 WITH INDEX 2*<I3*N(l)*N(2) + I2*Nd>+m + l AND THE IMAGINARY PART 

— r — i-N — iij. — t r.r.AT ion t mmf d t a tf i y foi iruiwc 

C IF THE FOUR IEP SFRIES IS REQUESTED, ARRAY A IS REPLACED BY 

C X( J1 ,32. J3 LsSUM A ( K 1 ,K2 . K3 ) *W1**(K1*J1 1 *W2**( K2 * J2 ) *W 3* »( K3* J 2 J .... 

C SUMMED OVER K2=0,N(2)-1, K3=C.N<3)-1 

- C— , — WHERE- W-UAUI.fcr.IH~ ROOT OF. UJNIXY. . 

-X M t- I -L. I =X «.2 «- 3 . WHF . RF Nt I \ =?T* Mi-I.l 1.S .THE _N.Qa. QE PISa IN.. THE I-TH. DIM. 

0 THF DIMENSION OF A IN THE CALLING PROGRAM SHOULD BE TWICE TEE 

0. NUMBER CF COMPLEX ELEMENTS OF .THE LARGEST A ARRAY TO BE PROCESS- 

0 FD 

_E THE XQMPXEX-Jt . » S. -ARF S T 0 R F D IN THE SAME MANNER AS A. 

0 

— C 1 E— I HE _£ £ UR..I E R TRANSFORM IS RE QIIF ST EQ, ...THE ARGUMENT A IS TAKEN 

C TO BE X AND IS REPLACED BY THE ARRAY A SATISFYING THE FOURIER 

C SERIFS* ... 

_C l FT MT=MAXf M(i) ,ML2),M(3) |-2, NT=2**MT, WITH M BEING THE M 

C GIVEN WFFN THE TABLES ARE SET. 


C SI J )»SIN( J*PI/(?*NT I U J * 1,2. 3,.. .NT-1. 


C I N V ( J+ 1 ) =WORD CONTAINING BITS OF J IN INVERTED ORDER IN ITS 

X RJXti3MU3 T^NI....BT IONSa. F.QR J = D ..t_1.2..e. .NT -1. ... 

C 

— C LEI T F S=0 TO SET UP S I N. AND I NV TABLES. 

C I F S = + 1 TO SET UP SIN AND I NV TABLES AND 00 FOURIER SERIES. 

X IFS =— I TO SFT LJg> S1N.ANQ I.N.V TABLES AND D.O.f OUR I ER^.TRANS FORM. 

C. I F S=+ 2 TO DO FOURIER SERIES ONLY. 

C IF S TO DO FOUR TER TRA NSF CR M O NLY. 


_C ONE DOFS NOT HAVE TO REPEAT THE CALL TO ’HARM* WITH IFS=Q,+1,-1 

C. IF ONF DOFS NOT CHANGE THE MAXIMUM Me 

C T FF RR =0 TF THF ARGUMENTS M ARE C. K« 

__£ 

f. T F FR P = 1 TF THERE IS AN FRROR IN CALLING 'HARM* 

C IF TFS=*. + 1 .-1 . TT MFAN S THAT THE MAXIMUM M IS GREATER THAN 20 

C. OR LESS THAN 3 

C IF JFS = + -2, IT MEANS THAT 4. $UFF IC IENTLY LARGE SIN AND INV TABLE 

C HAS NOT BEEN COMPUTED,. ONE MUST CALL 'HARM* WITH IFS=G,+-1 AND 

£ WITH A MAX M < I ) GREATER T HAN OR EQUAL TO THE MAX M( I ) FOR WHICH A 

C FOURIER TRANSFORM IS TO BE COMPUTED. 

£ 

C I FFR P =— 1 IF ONE IS CALLING CN 'HARM* WITH IFS=0,+-1 TO COMPUTE 
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f 


5 IN . T N V TABLES _ WHI CH IT .A L READY 

HAS 


c 

r. 


CALL TO HARM WITH THE SAME MAXIMUM M 



10 

i? 

IF ( TABS! IFS1-1 1 900,900.12 
M.TT = MA XO (Mil 1 .M( 21 ,M( 3) I -2 



R00T2 = SORT ( 2- 1 

IF ( MTT-MT > 14.14.13 


n 

_i 

I FER R = 1 
return 




14 

IFFRP=0 
Ml =Mfl 1 





M2=W(2 1 

V3=M( 31 . 





M1=?«*MI 

N2=2**M2 





N3=2**M3 

TF (TFS) 16,1,20 



r. 

16 

TO CAL CUI..ATF TRANSFORM REPLACE A 

NTOT ■ = __ .... 

BY 

CCNJGIAl/N 



FN = NTOT 
00 18 1=1. NTOT 




.18 

A(2*I-1) = A ( 2* I - 1 ) /FN 
A( 2*1 1 = -&(?*!) /FN 




?r> 

NP ( 1 1 =N1*2 
NP( 2 1 = NP < 1 ) *N2 





NP I 3 )=NP( 21 *N3 
HO 250 10=1.3 





TL = NP ( 3 1 -NP ( ID) 
TL1 = IL+1 





MI = M( 10 1 

IF (MI » 250.250,30 




80 

idif=np< mi 

KR I T =NP ( TO 1 





ME V = ?*(MT/2) 

IF (MI - MFV 16C.6C.40 



r. 

_ 40_ 

m is non- no 1=1 case 

KRIT=KRT T/2 





Kt =K RI T-2 

00 50 1=1 .TLl.IOIF 





KL A ST=KL+ I 

nn 50 K=T » KL A ST ,2 



_j i 


KD=K4KRIT 

DO ONF STEP WITH L=1,J=0 



c 

r. 


A ( K 1 = A ( K 1 + A ( KO 1 
A ( KD »=A( K )-A( KOI 



r. 


T = A ( KO 1 





A ( KD ) = A ( K 1 — T 
_ A(K1=A(K1+T 





T= A ( KD+1 1 
A ( K 0 4 1 ) = A ( K4 1 ) -T 




50 

A(K4] >=A(K4l ) 4 T 
IF (MT - 11250,250.52 



c 

5? 

(FIRST =3 

OFF - JLAST = 2**(L-2> -1 





.11 A S T = 1 





24 








M TS FVFN 

40-im«-a-4 

ji A^T=n 

70 on 240 L=l.FIRST,MI 
J JO I F=KR T T 


KL=KRTT-2 

— &a-cm~4*a 

on rc i*i«rn«iniF 

— JCLARXs T +KL — 

nn flr k=i,kla<;t,2 



K2=K 1 + KB IT 

K3=K 2-4-KB IT 



nn TWO STF P S W ! TH J= 0 



A(K t=A (K)+A(K 2) 



A ( K 1 ) = A ( K 1 ) +A ( K 3 ) 
_A-CJC34^UCU -ACK3J 


A< tLL=A. l tCUJU K.U— 
A ( K 1 »=A ( K )-A(Kl » 


A(K3)=A(K2)-A{K3>*I 
T=A ( K2 ) 

AUC2X=JUJa-T 

A(K l=A(K )*T 

-laAXK-ZUJ 

ACK2+1 ) =A ( K+l >-T 
ACK+1J = A ( K+ LULL. 


_ T.-ALK31 

A ( K 3 ) = A ( K 1 I — T 

A.(.LlL=AjLK.LLtX. - 

T = A ( K 3+ 1 ) 

AUO+L) -A (Kl + l )-T 

A(Kl+l »=A(Kl+l 1+T 


T= A ( K1 ) 


A ( K ) =A ( K ) + T 

T = A(K1 + 1) 

A ( K 1 + 1 ) = A ( K + 1 ) — T 
A ( K + 1 ) = A ( K + 1 L+T 


R=.-A.LO±ll._ _ 

T = A ( K 3 I 

A 1 K 3 1 = A i K 2J - R 

A(K?»=A(K2)+R 

..... A(K3+n=A(K2 + l J -I 

RO A ( K 2+ 1 l = A ( K 2+ 1 )+T 

IE LJI ASJJL _23JLjlZ25j,.82_. 

R2 J J = J Jf< IF +1 
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o <_ cj 


c 

nn fcr j=i 

II AST = IL +JJ 

• 


nn 8 5 I = JJ. I LAST. ID IF 
K1.AST = KL+I 

nn PS K = I .KLAST .2 

K1 = K+KRIT 

K ? = Kl+KRIT 
J <A = K ?+KR I T 

C 

.... C 

l FTTTNG V» = ( 1 + 1 ) / ROOT2 »W3 = ( — 1 + 1 1 /ROOT2.W2=I . 
.AIK.) =A(K)+A<K2 ) * I 



C 

C 

A(K2 )=A( K)-A(K2 >*I 

A (K1.L=A ( KL.) *h+A( K3) *W3 . _ 



C 

..... -C- 

A(K3)=A(K1)*W-A(K3)*W3 



c 

. c 

A(KI=A(K)+A(K1 ) 
_ AIK1)=A(K)-A(K1) 



c 

L 

A(K2 l=A(K2)+A(K3)*T 

. A(K3)=AIK2)-A(K3 ) * I ... 



n 

. R =“A (K2+1L. 




T = A ( K ? ) 

AIK?) = A ( K l-R 




A ( K ) = A ( K ) +R 

A( K2+1 ) =A ( K+l l ! - T 



c 

A( K+l * = A 1 K + l ) + T 




A WR = A ( K l )-A(Kl+l) 

A W I = A(K1+1) + A(K1 ) 




R=— A (K3 )— A( K3+1 ) 
T=A( K3I— A(K3+l ) 




A(K3 ) = (AViR— R) / Rf)nT2 
A(K3+1)=(AWI-T >/ROnT2 




A ( K 1 )=( AWR+RI/RnnT? 
A(K1+l)=(AWl+T)/RnnT2 



T = A ( K 1 ) 

A ( K 1 )=A( K)-T 

A( K >=A( K >+T 
T = A ( K I + 1 ) 


A(Kl+l)=A(K+l)-T 
A( K + l )=A (K+l ) + T 




R=-A(K3+I ) 
T = A C K3 ) 




A ( K 3 ) =A ( K? )-R 
... .AC K ? ) = A ( K2 ) +R 




A ( K 3+ 1 ) = A I K 2+ 1 ) — T 
85 A C K ? + l ) =A ( K 2+ 1 ) + T 




I F C .11 AST-1) 235«?35«90 
90 JJ= JJ + JJOIF 




c. now nn the remaining jmj; 
nn ?30 ,i=? . jl a s t 


FETCH W • S 

DE F- W=W**INV(J), W2=W«*2. W3=W**3 

96 I=TN\/’f.l+l) 
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• — <5R. IC=NT~l 

W( 1 ) =S ( If I 

w( 2 > =s< 1 1 

T2=2*I 
I ?C = NT— -1-2 

if u 20 ) 120 , 110 . ion 
C 2*1 IS IN FIRST QUADRANT 

ion w?n » = s( i.?r. » 

W2 ( 2 ) = S ( 12) 

GO TO 1 30 

110 

W2X24.=1-- 

GO TO 130 

C - - 

r 2*1 IS IN SECOND QUADRANT 

120 I 2C.C - = 1 20 NT. 

T 2C a— T 20 

W.2X l_Ls.-SXL2.ILl 

W 2 f 2 ) = S ( 1200 ) 

130 13 = 1 + 12 

T 3f. =NT- T 3 

- I F ( I 30 1 160.150.140 

C 

~C. - L3-..IA E I£.S. L-- QU ADlP.AXJ. 

140 W3( 1 l=S( T 30 ) 

W312 )=S< 131 

GO TO 200 

150 W3C 1 X=0* 

W3< 21=1* 

GO-in.. 200 

0. 

160 I 30. 0 = 130 4.NT 

I F ( I 300 I 190.180,170 

C. _ 

0 13 IN SFOOND QUADRANT 

1X0,3 30 =-130 

W3< 1 )=-S( 130) 

— „ W3.J.2-) = S ( 1 200.1 

on TO 200 

1 80 -W3LLll.=-l a _ 

W3( 2 ) =T, 

GJCL.TJC._ZOC 

0 

--C-- - _ 3*1 I N THIRD .QUADRA.MT 

190 I 3 0 0. 0 =N T + 1 3 CC 

— 130.0 = - T 300 _ __ 

W3( 1 ) =— SCI 3000 ) 

M3±21a^S±L2LU 

200 TLAST=IL+JJ 

— DO.. 2 20 J = J J.t.I. LA ST , I D I F ..... ... 

KLA ST = KL4 I 

DO 220 K = I *_KL A S T,2 

Kl =K4KBI T 

L 2 =ja±joa,Li 

K 3=K 24KB I T 
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P n\n 


c nn two steps with j not o 

C A ( KJ ~A.( K ) + A ( K2 1 *W2 

C AC K2 * = A I K » — A C K 2 * *W2 

C A ( K 1 > =A C K 1. ) *W+ A ( K 3 L*M3_. . _ 

f. ACK31=A(K1 )*W-ACK3t *W3 

f. A C K ) =A ( K ) +A C K 1 I 

C A (K 1 ) =A IXL-AiKli 

C AC K? )=A(K2 )+ACK31*I 

C A C K 3 . 1 =A f K2.1-MK.3 ) * I _. ... _ 

r. 

R.=ATXZ)*I&2 U1.-M.K2+ li*J*2{.21 _ 

T=A(K2)*W2C2»+A(K2+1)*W2(1) 

A(K2)=A(K)-P 

ACK)=A(K1+R 

ACK241 )=A(K+1)-T 

ACK + 1 1 = A ( K + 1 1 + T 

R = A ( K 3 )*W3 ( 1 )-A( K3+1 >~*W3( 2 > 

T = A ( K 3 ) *W3 ( 2 ). + A (K3 + 1 ) #.W3 ( 1 ) 

AWR=A<KU*U(l)-A(Kl + l)*W(2> 

AWT = A ( K1 )*WC 2 > + A CKl + l ).*W( l ) 

A < K3 1=AWR-R 

. A..( K.3 + 1J =A W.I-I. .... .... _ „ 

AC K1 ) =A WR+R 

A C K 1 + 1 J = A W I + T 

T = A C K 1 ) 

AC K 1 ) = A ( K 1 — T _ 

ACK »=A(K l+T 

T = A ( Kl+1 > _ _ 

ACK 1+1 j = A ( K + 1 )-T 

ACK + 1 ) =ACK + 1 ) +T 

R=-A(K3+1) 

T = A C K 3 1 

A(K3)=A(K21-R 

A(K21=A CK2)+ R 

ACK3+1 1 = A C K 2 + 1 )-T 

. 220 A C.K2+J 1 =A ( K 2 + 1 ) + T 

0 FNO OF I AND K LOOPS 

? 3Q J J -J JD I F + .I J 

0 END OF J-LOOP 

235 J LAST=4 *Jl AST+3 

240 CONTINUF 

C END OF L LOOP .. 

250 CONTINUE 

FND CF ID LOO P 

WF NOW HAVE THF COMPLEX FOURIER SUMS BUT THEIR ADDRESSES ARE 
C RIT-RFVFRSED. THE FOLLOWING ROUTINE PUTS THEM IN ORDER 

NTSO=NT*NT 

M3MT =M3-VT 

3 SO IF CM3MT1 370,360 ,360_ 

C M3 GR« OP F 0, MT 

3 60 I G O 3 = 1 

N3VNT=N3/NT 
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MINN 3= NT - 

GO TP 380 

C M3 L.FSS THAN MT 

3 70 T GO 3 =2 
N3VMT=1 
NTVN3=NT /N3 

...... -MI-WW3=N3 

380 J JO 3 = NTS0/N3 

M?MT=M?-MT s 

450 IF IM2MT147C.46C.460 

C M2 G R« OR F0« MT - - - - 

460 T GP? = T 

U2M HJjsU2JHZ 

MINN?=NT 

GO TO 480 - _ — - — - - - 

C M? LFSS THAN MT 

470 I GO? =... 7 - - - 

N 2 V M T = l 

NXVM?= MI/AL2... .... 

MT NN2=N? 

480 J.W?=NTSC/N? 

M I M T =M T - MT 

550 IF( MIMT) 570,560.560 . .... . - 

C Ml G R« OR fOr. MT 

560 tGOAs.1 : 

N 1 VN T=N 1 /NT 

- - MT NN 1 =.NT 

GO Tp 580 

0 Ml l.FSS THAN MT 

570 I GO 1 = 2 

... N.1A/N.I=..L 

NTVNl =NT/NI 

MT NN1 = N1 - 

580 Jjni=NTSC/Nl 

600 J .13-1 _ - 

J = 1 

HQ 8fi.0._ ^Lg.P-3 =JL jlN 3^/lhLI 

I PP3=I NV ( J J 3 l 

no 870 JP3-1.M I N.N3 

GO TP (610* 620 1 « IG03 

610 !P3 = INV ( J°3)*N3V.MJ 

GO TO 630 

A20..ae3=imt.LJi?3iy.MIVli3 

630 T 3= ( TPP3+TP3)*N2 

... 700 JJ?=1 

no 870 JPP?=1 « N2VNT 

IR E 2 = I N V I J JL2 J± 1 3 

no 860 JP2=1»MI NN2 

G O TO J 7 1C.i 7 2Q.lAjfi.Q2 .- 

710 T P 2 = I N V ( JP2 ) *N2 VNT 

GO TO 7.3.0 

720 TP?=INV( JP21/NTVN2 

7 30 12 = ( TPP2+ I R.2J * N 1 

800 J.M =1 

D Q . 8 6 0 „ JRR 1 = L,Ji mil 

T P P 1 = 1 N V ( J J 1 )+ T 2 
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r> ,n jn poh n 


DO 8 5 C. . JLR J = 1 ». M I N.N 1 .. 

no to ( pn ,820 > , igoi 

810 IP1=INV( JP1)*N1VNT 
GO TP 83 r 

820 IP1=INV< JPD /NTVN1 
830 !=?*( IPP1+IP1I+1 

_ IF. ( J- I I . 841.88 5,845 

840 T=A( I > 

Ai.Il.sA.Lil 

A ( J ) =T 

T = A (1+ 1 ) 

A ( I + l » = A C J+l 1 

A.LJ + .1 L=T. 

845 f.ONTINUF 
850 J =J+2 
860 JJ1=JJ1+JJD1 
C FND OF JPP1 AND JP2 

870 JJ?= JJ2+ JJ02 

£ .FND. _QF, J PP7..AN.D J P.3 LOOP S 

880 J J 3 = JJ? + JJD3 
C END OF JPP3 LOOP 

IF(IFS) 882,1,1 

0. DOING TRANSFORM,, P F PLACE A BY CONJG(A). 

882... 00. 8 84 .1=1., NIPT. .. . ... 

884 A( 2*1) = - A ( 2* I ) 

GO TO 1 ... 

R F TURN 

THF FOLLOWING PROGRAM COMPUTES THF SIN AND INV TABLES 

90 0 M T^'f^TcTF(Ti“ M ( 2 l", M"nT V ~2 

M T = M A XO (2 » M T ) 

90 A- IF ( MT-20 1906,806, *505 
005 I F F RP = J _ 

GO TO 1 

R FTURN ___ 

906 YfFRR=0 

N T = ? **MT 

N T V 2 =M T / 2 

SFT__ lfP_ STN T ABl F 

THFTA=PIF/?**(L+1> FOR L=1 

9 10 THFTA = „ 78 53581 634 _ 

C J S TF P= ? * *( MT— I + lT F~ OR L = 1 

JSTFP_=NT 

0. JD IF=?**( MT-L > FOR L=1 

JD T F =NTV 2 
S< JOIF )=SIN(THFTA> 

DO 950 L = ? » MT _ „ 

THF T A = TH FT A 7 ?7 

J S.TF P?= JSTF P 
JSTFP=JDIF 
JDI F=JSTFP/2 
S( JDIF)=SIN(THFTA) 

ifj -NT- JD I F 

S ( JO. 1 ) =CT S~( THF T A ) 
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- -JL-A JFJ?-2 - 

TFCJLAST - JSTFP) 950,920,920 
920 nn 940 J=JSTFP,JLA5T, JSTFP 
.IC = NT-J 

J0 = J+JDIF - . .. 

940 jm»s( j)*s< jr.n+s(jDiFi*s(jci 

9*0 CONTINUE- - 

SFT UP I NV ( J ) TABLE 
96 C Mil F yp = NTV? 

0 MTLFXP=2**< MT-L »o FOR L = 1 

IMiFXPsI 

C \ Ml Fyp = ?**(l-1 )« FOR L= 1 

TNVm=0 
on 9 BO L=1,MT 
TNVUMlFXP-rlJ = Mil EXP 
00 97f J=2,LM1FXP 

JJ = J*LM JFXP 

970 IMV( J J I =INV( J) 4MTLFXP 

MTJ F XP = MTLE XP/2 - 

980 i Mipyp^Mieyp*? 

!F( IFS) 12.1,12 

0 RFTURM 

FND ... 
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